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We introduce an exactly solvable model for glassy dynamics with many relaxational modes, each 
one characterized by a different relaxational time-scale. Analytical solution of the aging dynamics 
at low temperatures shows that a nonequilibrium or effective temperature can be associated to each 
time-scale or mode. The spectrum of effective temperatures shows two regions that are separated by 
an age dependent boundary threshold. Region 1 is characterized by partially equilibrated modes that 
relax faster than the modes at the threshold boundary. Thermal fluctuations and time-correlations 
for modes in region I show that those modes are in mutual thermal equilibrium at a unique age- 
dependent effective temperature 0(s). In contrast, modes with relaxational timescales longer than 
that of modes at the threshold (region II) show diffusive properties and do not share the common 
temperature 0(s). The shift of the threshold toward lower energy modes as the system ages, and 
the progressive shrinking of region II, determines how the full spectrum of modes equilibrates. 
As is usually done in experiments, we have defined a frequency-dependent effective temperature 
and we have found that all modes in region I are mutually equilibrated at the temperature 0(s) 
independently of the probing frequency. The present model aims to explain transport anomalies 
observed in supercooled liquids in terms of a collection of structurally disordered and cooperative 
rearranging mesoscopic regions. 

I. INTRODUCTION 

Finding a general theory for nonequilibrium systems constitutes one of the major challenges in modern physics. 
Despite of the many different existing approaches to nonequilibrium phenomena 1 ' 2 , at present there is not yet a 
general theory applicable to systems arbitrarily far from equilibrium. During many decades, physicists have studied 
different sorts of nonequilibrium systems. A common strategy has been to adapt well-known thermodynamic concepts 
to the nonequilibrium regime with the hope that the most fundamental assumptions still hold in a varied range 
of nonequilibrium conditions. One of the concepts that physicists have tried to rescue repeatedly is that of local 
equilibrium 3 and the possibility of defining a nonequilibrium temperature 4 . Different definitions of a nonequilibrium 
temperature have been explored for nonequilibrium liquids 5 , steady-state driven systems 6,7 , granular media 8 and 
glassy systems 9,10 . Admittedly, the concept of a nonequilibrium temperature in realistic systems still remains elusive, 
mainly because of the scarce number and difficulty of the experiments that have to address fundamental predictions 
by the theory 11,12 . 

Compelling results have been obtained over the past two decades in the context of glassy systems. Several concepts 
(e.g. mutual equilibration) have been rationalized by borrowing mathematical tools used in the context of spin glasses 
(such as reparametrization-time invariance or supersymmetry 13 ). Glassy systems are non-stationary systems which 
relax very slowly toward equilibrium 14 . Relaxational processes can be characterized by the age of the glass, which 
is the time elapsed since the system was prepared in the non-equilibrium state. A key signature of aging is the fact 
, that correlations and responses decay in a time-scale which is roughly proportional to the age of the system, thereby 
showing strong and complicated hysteretic effects 15 . Theoretical studies in mean-field spin-glasses have shown 9 that a 
nonequilibrium temperature can be defined in the aging regime. This entails an extension of the fluctuation-dissipation 
theorem (FDT) to the aging regime. This nonequilibrium temperature (usually called effective temperature in the 
glass jargon) is defined by, 

/ 9C(t,s) \ 

k B T cS (t, S) = [ 1 t > S , (1) 

where C(t, s) is a generic two-times correlation function and G(t, s) is the corresponding response of the system to an 
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external perturbation applied at a given previous time s. At first glance, (1) is just a useful way to quantify deviations 
from thermal equilibrium (the equilibrium regime is characterized by T c g(t, s) — T). 

In experiments 11 the effective temperature is usually defined in terms of quantities measured in frequency space 16 : 

where is the probing frequency, S(Cl, s) is the power spectrum of fluctuations at time s and x"(^> s ) is the dissipative 
part of the response. In equilibrium T c g = T independently of the probed frequency fi. 

Over the past years the validity of this parameter from a statistical and thermometric point of view has been widely 
studied 17,18 . Yet, there are still several debated issues. Of the utmost importance is to understand under which 
conditions the effective temperature governs the heat flow between different aging systems, whether a zero-th law 
is applicable between different systems sharing a common effective temperature 19 and whether such temperature is 
experimentally measurable 12 . More speculative but not less important is whether we could use such measurements to 
trace back the age of a man-made glass, i.e. to approximately determine the time elapsed since it was quenched. Our 
capability to answer such a question might help us to understand the limitations and difficulties in the experimental 
measurement of effective temperatures. 

No doubt that many aspects related to the noncquilibrium properties of glasses are tightly dependent on the 
microscopic details of the system under study (a polymer blend, a structural glass former, a disordered magnet or 
a dirty superconductor). However, generic and universal aspects of the noncquilibrium thermal properties in glassy 
systems are expected beyond the specific properties of the system under study. This has motivated the appearance 
over the past decades of several phenomenological approaches to the glassy state, from the old free volume theories to 
the Adam-Gibbs-DiMarzio entropic theories 20 , the free-energy landscape approaches 21 until the more recent mosaic 
pictures 22,23 . 

A major question that has recently attracted much attention in experimental studies of supercooled liquids is the 
microscopic origin of the non-exponential behavior observed in the relaxation functions. According to some views, 
a glass is spatially homogeneous and the relaxation of local regions in the system is intrinsically non-exponential. 
Another view supports the idea of a glass as an intrinsically heterogeneous system made out of a collection of 
cooperative mesoscopic regions in the system, each one showing exponential relaxation with its own characteristic 
relaxation time. Ergodicity also implies that regions must have a spectrum of finite lifetimes and that these regions 
intermittently switch from mobile to frozen. The superposition of the relaxation of the many local regions gives rise 
to the non-exponential behavior observed in the relaxation of bulk quantities. Concomitant to the experiments there 
has been an upsurge in the theoretical study of heterogeneities in simple models using analytical methods 24-26 and 
numerical simulations 27 . In general these models support the heterogeneous scenario in agreement with what has 
been observed in experiments 28-30 . If the heterogeneous scenario is valid then one can ask whether it is possible 
to define a mesoscopic effective temperature locally defined for each region and intermittently fluctuating at time 
intervals roughly equal to the lifetime of the region. In such case it seems more reasonable to talk about a fluctuating 
effective temperature field rather than a single well defined value. Recent AFM 29 , fluorescence spectroscopy 30 and 
Nyquist noise measurements 11 ' 31 seem to endorse the heterogeneous scenario. Moreover, Nyquist noise experiments 
purport evidence of frequency dependent effective temperatures (2) which can be interpreted in terms of degrees of 
freedom that have different effective temperatures and not being in mutual thermal equilibrium. 

Inspired by these facts, we introduce an exactly solvable model with glassy dynamics, the Disordered Oscillator 
Model (DOM). The model is an extension of a previous one introduced in Ref. 32 as a simple model for a glass with 
slow dynamics in the presence of entropy (rather than energy) barriers, the so-called Oscillator Model (OM). The OM 
represents a simple example of cooperative slow dynamics due to entropy barriers where kinetics is constrained by 
the dynamical rules rather than by the Hamiltonian. In the present paper we show how a mode dependent effective- 
temperature can be defined in the presence of many different time-scales. This is naturally accomplished by adding 
structural disorder in the OM (a similar strategy has been followed for the Backgammon model 33 ). In the presence 
of disorder the system preserves aging dynamics at low temperatures and a mode-dependent effective temperature 
which is higher than the bath temperature drives the relaxation of the system toward equilibrium. The scenario 
of relaxation in this model is as follows: relaxational modes can be split in two regions (I and II) separated by a 
threshold boundary. Region I contains mutually thermalized modes at an age-dependent (denoted by s) but frequency 
independent effective temperature 0(s) that lies above the bath temperature T. In contrast, modes in region II are 
non-thcrmalizcd and diffusive. Relaxation takes place by the progressive mutual thcrmalization of degrees of freedom 
in region I to the single temperature 0(s). As the system ages, the threshold region shifts, O(s) decreases further 
and more modes enter in region I. Full thcrmalization occurs when O has relaxed to T and region II has vanished. 

The paper is organized as follows. In Section II we introduce the model. In Section III we present the formal 
solution of the model. In Section IV we derive the dynamical evolution for the energy, correlations and responses. 
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In Section IV we also analyze the behavior of the non-equilibrium effective temperatures for the different modes in 
the system. In Section V we derive the effective temperature using fluctuation relations. Section VI summarizes the 
results and Section VII presents a discussion of the implications of this work to our understanding of the glassy state. 
Technical aspects are described in the Appendices. 

II. THE DISORDERED HARMONIC OSCILLATOR MODEL (DOM) 




FIG. 1. Schematic picture of a system of independent oscillators with different springs constants. The thickness of the 
springs represents different values for the spring constants. The motion of the oscillators is not independent due to the coupling 
induced by the Monte Carlo dynamics (represented by the top rectangular box with the heading "dynamics"). 

The original OM model was introduced in 32 as a simple entropic model for a glass. It consists of N identical and 
uncoupled harmonic oscillators with Hamiltonian: 

i=l 

where K is the stiffness or spring constant and the Xi are the positions of the oscillators. The dynamics of the model 
is stochastic of Monte Carlo type. A Monte Carlo trial move is defined as follows. The position of all oscillators 
are simultaneously shifted by Xi+ri/ y/N where r» are random variables extracted from a Gaussian distribution with 
zero mean and variance A 2 . The global move is accepted according to a transition probability W(AE) which satisfies 
detailed balance: W(AE) = W(-AE)exp(-/3AE), where AE is the change in the energy (3) and (3 = \/{k B T) 
(from now on we set the Boltzmann constant equal to 1, ks = 1). We use the Metropolis algorithm for the dynamics. 

W(AE)=mm[l,exp(-(3AE)} . (4) 

Note that whenever a move is accepted, all oscillators arc updated in parallel which means that the motion of all 
oscillators is cooperative and coupled by the dynamical rule. Therefore, although they are non-interacting in the 
Hamiltonian (3), the dynamical rule couples all the oscillators in a non-trivial way. 

At low temperatures the OM shows typical non-equilibrium features of glassy systems such as aging in the correlation 
and response functions. Indeed, at very low temperatures the acceptance of the Monte Carlo updates decreases very 
fast (like 1/t). In this way relaxation generates entropic barriers which make dynamics extremely slow (the energy 
decreases as 1/logt). The simplicity of this model, however, makes it exactly solvable. For a compilation of other 
results see Refs.( 16 ' 34 ' 35 ). The OM has also been extended to include non-harmonic generalized potentials 36 . 
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In the present work we add different time-scales to the system by introducing disorder in the different degrees 
of freedom that enter into the Hamiltonian. The disordered oscillator model (DOM) is defined by the following 
Hamiltonian: 

1 N 

H Q = -KJ2^ • (5) 

i=l 

where Kti is the spring constant of each oscillator. As in the original model, the system consists of a set of non- 
interacting harmonic oscillators evolving in parallel according to a Monte Carlo dynamics (4). Again, the dynamics 
is cooperative and couples the evolution of all oscillators. In this disordered version, the spring constants of the 
oscillators have a probability density given by, 

■9( £ )-^E^- £ ) • ( fi ) 

i 

We assume that et > 0, Vi. The OM (3) is a particular case of the DOM with g(e) = 5(e — 1). In Fig. 1 we show 
a schematic picture of the model where a set of different oscillators move coupled according to the dynamical rule. 
The inclusion of disorder in the DOM allows us to explicitly consider nonequivalent degrees of freedom in the system 
having different relaxational properties. 

Thermodynamics of the DOM is trivial. The free energy per oscillator / is given by, 

1 9-7T 1 f°° 

/ = ~2 l0g W + 2 I ^( e ) lo §( £ ) ' ( ? ) 

where g(e) is defined in (6). This expression is valid for any finite N if q > 0,Vi. In the large N limit, g(e) may 
become a continuous function of e. The last integral term in (7) can diverge depending on the behavior of g(e) when 
e — > 0. Whatever value this last term takes, it only represents an additive constant to the total entropy. It may be 
taken as defining the entropy of a reference state without influencing the rest of thermodynamic quantities. As we will 
see later, the behavior of the oscillators with coupling constants e — > determines the global relaxational properties 
of the system. 



III. FORMAL SOLUTION OF THE MODEL 



In this section we derive the dynamical equations for the energy and the acceptance rate of the system. To this 
end we must compute the probability distribution of having a change SE in the energy E — Hq/N (5) in one Monte 
Carlo step. In a Monte Carlo step the N oscillators are moved according to the rule: 

Xi — ► Xi + -j= Mi , (8) 
where rj arc random variables Gaussian distributed with zero mean and variance A 2 . The change in energy is given 

by, 

N 

(9) 



/ 2 \ 



In what follows we identify ^ j=1 by J2 i without risk of confusion. The probability of having a change SE can be 
written as, 



p ^=Lj( se - k ^-{^+0 (n 



dr, 

■. exp 



(10) 

Using the integral representation of the Dirac's delta function and performing the Gaussian integrals we obtain: 

1 / (SE — ka2 £ ) 2 \ 

P(SE) = - exp — - — Trrp. — , (11) 

1 ' V4^Q 2 A 2 I 4^g 2 A 2 I v ' 
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where we have defined the following quantities: 



^=<^E e ^> ■ ?= ^E e - ( 12 ) 

i i 

where the brackets < . . > denote an average over initial conditions and dynamical histories of the system. 

Note that (11) is analogous to that obtained in the OM 32 (where e 4 = 1 ,Vi, and Q k = KE,Vk). In the DOM the 
energy is given by E = Q\. Comparing to the OM (11) now depends on Q 2 rather than on the energy E. As in the 
case without disorder, the evolution for the energy is given by the expression: 

dE f 00 

— = / dxxP(x) + / dxxP(x) exp(-/3x) , (13) 

&t J-oo JO 

where we have used (4). Straightforward calculations give, 

^ = T l er fc(a) + l 1 ) ^ {t) ) ' (14) 

where we have defined the following quantities 

a c = -KA 2 e, a {t) = —======, eric{x) = -= / exp(-y 2 )dy, (15) 

2 ^AKA 2 Q 2 (t) Vtt J x 



f Q2 (t) = exp (-(3a c (l ^y)) crfc (a(t) 



( /3a c 



(16) 



It is easy to see that for e, = 1, Vi, the results of Ref. 32 are recovered. The acceptance or average fraction of accepted 
moves is 

/*0 />oo 

A(t)= dxP(x)+ dxP(x)cxp(-(3x) . (17) 



Inserting (11) in (17) we get: 



A(t) = l(erfc(a(t))+/ Qa (t)) . (18) 

To solve (14) and (18) we note that the equations for the energy and the acceptance are not closed because they 
depend on the quantity Q2{t). The equation for Q 2 turns out to depend on Q 3 and Q 2 - In general, the equation for 
Qk depends on Qk+i and Q 2 . In Appendix A we present the dynamical equations for Q k . The final result (A7) is: 



dQ k _(Ktfe k a c Q k+1 \ m+ (Q^\dE 



dt V 2 Qi ) \ Qi ) 9t 

where we have defined the moments of the coupling distribution: 



1 f°° 

. I " k 9(^ ■ ( 2 0) 



In what follows we will identify e = t\. We must notice that the hierarchy (19) can be solved in the same way as in 
the case of the spherical Sherrington-Kirkpatrick model 37 or the Backgammon model 38 by introducing a generating 
function and solving the resulting partial differential equation by the method of the characteristics. We will not pursue 
this strategy here, rather we will analyze the simpler adiabatic approximation. 
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A. Adiabatic approximation: the zero-temperature limit 

Equation (14) can be solved at zero temperature by means of an adiabatic approximation. It states that some 
quantities evolve much slower than the energy so we can consider them as partially equilibrated over the constant 
energy surface. At equilibrium the derivative of the energy (14) is zero and Q 2 q is given by: 

Q° 2 q = |?=^e . (21) 

The adiabatic approximation becomes asymptotically valid for long enough times. It relates Q2(t) with the energy 
E(t) by assuming that the system has partially equilibrated over the hypersurface of energy E(t): 

Q 2 (t) = E(t)e (22) 

Within the adiabatic approximation (14) is closed and can be solved analytically. In fact, this equation is identical 
to the equation of the energy in the OM model. As was done in Ref. 32 the time-evolution equation for a(t) (15) at 
zero temperature is: 



da _ cxp(-q 2 ) 

where we have expanded the error function in the limit a. — > oo. For large times, the parameter a grows logarithmically 
as a function of time: 

f It 1 It \ ^ 

a(t) ~ (log(^) + ^ log(Iog(-^))J , (24) 

therefore, the acceptance is independent of the disorder distribution g(e) and decays as, 

1 

4*log(-%) 



Mt) * 2t , , (25) 



plus subdominant logarithmic corrections. The energy decays logarithmically in time and depends only on the mean 
value of the distribution, c. 



KA 2 e 1 

{t) " 16 Iog(^) + ilog(Iog(^)) • (26) 

In Fig. 2 we show the time evolution (measured in Monte Carlo steps) of the energy for a Gaussian distribution of 
the disorder. The dashed line filled with circles is the asymptotic result (26) which is in agreement with the numerical 
solution of (14) (continuous line) for asymptotically large times. In the inset of Fig. 2 we show the acceptance. Again, 
circles are the asymptotic solution (25) while the continuous line is the numerical solution of the exact equation (18). 
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FIG. 2. Zero-temperature time evolution of the energy for a Gaussian distribution with mean ? = A = 5. The time is 
measured in Monte Carlo steps and the energy is dimensionless. The dashed line with circles is the asymptotic solution (26) 
while the continuous line is the numerical solution of (14). In the inset we show the acceptance. The dashed line with circles 
correspond to the adiabatic solution (25) and the continuous line is the exact solution of the numerical equations. To solve 
(14,18) we have integrated up to 500 terms in the hierarchy (19). 



From (21) and the adiabatic approximation (22) we can define a time-dependent effective temperature (that we 
denote as 0(f)) given by: 

Q(t)= 2Qg)_ 
e 

At this stage of the calculations this is just a definition which denotes the temperature that we could assign to the 
system if it were partially in equilibrium over the hypersurface of constant energy E(t) (meaning that it reaches 
the maximum entropy state compatible with that energy). In the following sections we will show how this effective 
temperature drives the relaxation of some degrees of freedom in the system. In addition, we will show that 0(t) can 
be exactly derived by assuming that energy fluctuations are described by a non-equilibrium measure reminiscent of 
the Edwards measure introduced in the context of granular media 39 . 

IV. ANALYTICAL SOLUTION OF MODE-DEPENDENT QUANTITIES 

In this section we analyze the relaxation of the different degrees of freedom in the system. We define the mode- 
energy density which is the energy density of those oscillators with a given spring constant value Ke, and compute 
its time evolution. Then we compute the correlations and responses for the different modes. Finally, we study the 
effective temperature for the different degrees of freedom in the system. 

A. MODE-ENERGY DENSITY RELAXATION 

A physical quantity which provides deeper insight into the behavior of the system is the mode-energy density w(e, t) 
defined as follows: 

^,*)=<^5>(ei -*)*?(*)> • ( 28 ) 
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This quantity corresponds to the mode-energy density of those oscillators with the same spring constant Ke. The 
total energy of the system can be written in terms of the mode-energy density as: 

/>OC 

E(t) = / ew{e,t)de . (29) 
Jo 

The characteristic relaxation time r(e) of a given mode is inversely proportional to Ke (see the computations in 
Appendix B). Therefore, w(e,t) is the energy density of those oscillators with a given characteristic relaxation time 
r(e) (B8). In what follows we consider the time evolution of w(e,t). The computation is outlined in Appendix A and 
the result (A13) is: 

dw(e,t) _ few(e,t)\ 8E + / KA 2 g(e) _ a c ew(e,t) \ ^ (3Q) 



dt V Q2(t) J dt \ 2 Q 2 (t) 

In equilibrium one-time quantities are stationary. So the derivatives of both the energy and w(e, t) are zero. From 
(14) and imposing stationarity for the energy, we obtain the relation (21) between the temperature and Q^ q - From 
(30) and imposing that the derivative of the w(e, t) is zero, we obtain the following relation: 

ew ^(e) = = Z g{e) , (31) 

valid for any e in equilibrium. Eq. (31) is the equipartition theorem for each energy mode e. It states that the average 
energy per degree of freedom is equal to . 

Next we consider the mode-density relaxation at T = 0. Eq. (30) can be solved in the asymptotic long-time limit 
using the results from the adiabatic approximation (25,26) by noting that ^ << A(t) and dW Ql'^ — ► yielding: 

ew(e,t) ~ 4^Q 2 (t) Ve , (32) 
which can also be written using the adiabatic approximation (22) as: 

ew{e,t) ~ g(e)E(t) -> < ^ {oT t ^ ^ . (33) 

This last result, which is valid for all values of e in the asymptotic (adiabatic) long time limit, is based on the fact 
that the energy is a monotonic decreasing function of time. On the other hand, for low enough values of e the first 
term in the rhs of (30) can be neglected and we get: 

d -^ = KA2 ^A(t)>0 e^O . (34) 

We conclude that, for long-enough times and for continuous distributions of g(e) with support on e-values extending 
down to e = 0, there must exist a value e* such that dw % — — 0, i.e. l9t "i^'*' > changes sign at e = e*. From (30,33,34) 
we can distinguish three different regions of e values depending on the sign of dW gl'^ '■ 

Region I : ew(e,t) > 9 -^Q 2 {t) ( ^—^ j — < ' ( 35 ) 

Region II : ew(e,t) < 4^Q 2 (t) ( \—be) — > ^fe^ > , (36) 



Threshold Boundary : e*w(e*, t) = ^V-Q 2 {t) ( iSe) — » 1 ' = ■ ( 37 ) 

Equation (37) defines the time-dependent threshold e*. Modes above the threshold are described by (35) and have 
a negative time-derivative of the mode-energy density w(e,t) while modes below the threshold are described by (36) 
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and have a positive time-derivative of the corresponding mode-energy density. For distributions g(e) whose support 
of e values does not extend down to zero it does not necessarily exist a value e* which satisfies (37). In such cases the 
two regions (35,36) may still exist but not the threshold boundary. 

At zero temperature the system is always out-of-equilibrium: the acceptance decays as A(t) w tl * t (25) and the 



energy decays as E(t) 
limit: 



log* 



(26), (these results being valid at leading order in 1/logf). Therefore in the long-time 



dE 



1 



A(t)a c dt " log* 



The asymptotic time-dependence of the threshold (37) is then given by the simple expression: 



(38) 



(39) 



From the dynamical equations (30) we can check that modes with e > e* verify the equipartition theorem at all times 
(32), 



ew(e,t) = ^-Q 2 (t) 



for 



e > e 



(40) 



so these modes are partially equilibrated at the temperature @(t) (27), i.e they have equilibrated over the hyper 
surface (in phase space) of constant energy E{t) = 

From (22), the value of Q2(t) decreases as the energy does. For distributions 17(e) with support on e- values extending 
down to e = 0, the threshold value e* decays asymptotically to zero as a function of time. The time decay depends 
in a non trivial way on the form of the distribution g(e) for e — > 0. For a distribution with finite weight at zero e the 
decay seems to be a power law (see below) . Unfortunately we have been unable to find a general analytical expression 
for such decay. 

By numerically solving (39) we can compute the time-evolution (in Monte Carlo steps) for the energy threshold. 
In Fig. 3 we show the time evolution of e* for a system that consists of 100 types of oscillators with g(e) — 1 and 
spring constants uniformly distributed in the range e = 0.01 — 1. It is clear that, as time goes on, the threshold energy 
diminishes and more and more oscillators partially equilibrate. As the system under consideration has 100 different 
types of oscillators then the solutions of (39) are discretized also. This is the reason why the time evolution for e* 
shows a stepwise shape. 




FIG. 3. Time evolution (in Monte Carlo steps) of the threshold energy value (39) for a system of oscillators with g(e) = 1, 
e £ [0, 1] for 100 spring constants uniformly distributed between 0.01 and 1. As time goes on, e* decreases meaning that modes 
with longer time-scales progressively equilibrate. 
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As the energy of the mode gives the characteristic time-scale of the corresponding oscillators (see Sec. IV B below 
and Appendix B, in particular (B8)), those oscillators which have a higher value of the spring constant Ke also have 
a shorter relaxational time-scale. As time goes on, the modes with e > e* are the first to partially equilibrate, while 
those modes with lower e equilibrate over longer time scales as they are cached up by the decreasing threshold value 
e* (39). 
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FIG. 4. (Color online) Zero-temperature time evolution (measured in Monte Carlo steps) of w(e, t) (30) for a system with a 
binary distribution of the disorder g(e) — (l/2)(<5(e — ei) + S(e — 62)) with ei = 1 and €2 = 0.1. We also plot w(e,t) coming 
from the corresponding equipartition relations (40) for the two modes of the system. The two continuous curves in the top 
correspond to ej while the two curves at the bottom correspond to £2 ■ The symbols correspond to the Monte Carlo simulation 
for a system with 5000 oscillators in each mode. The continuous lines correspond to the numerical solution of (30) and (40) 
corresponding to the exact and the adiabatic solution respectively. Different regions I (35) and II (36) according to the sign of 
the time derivative of w(e, t) are also shown. 

A summary of this behavior is shown in Fig. 4 for the case of a binary distribution of the disorder g(e) = (l/2)(S(e — 
e%) + 5(e — £2)). For this particular example we note that equation (39) is only defined for e = 61,62 so it does not 
have a well defined threshold solution e*. However, as we have already mentioned earlier, we can distinguish the 
two regions (35) and (36) corresponding to different signs of dW gl^ ■ This is shown in Fig. 4 where we show the time 
evolution of w(e, t) at zero temperature (30). We also plot the solution for w(e, t) using the equipartition theorem (40) 
for each of the two modes. We can clearly see that the oscillators with the largest value of e partially equilibrate very 
fast to the temperature 8(i) (27) (the two curves at the bottom), while the two upper curves deviate from each other 
even at long times. In Fig. 4 we also plot the results obtained from a Monte Carlo simulation (symbols) and compare 
them with the numerical solution of the exact equation (30) and the adiabatic approximation (40) (both in lines). 
The adiabatic solution reproduces well the exact solution at large values of e (region I) whereas it deviates from it at 
the lower values of e (region II) . The collapse between Monte Carlo simulations (symbols) and the numerical solutions 
(lines) is excellent. 

To sum up, in the non-equilibrium regime we have found two different classes of oscillators in the system: on the 
one hand there are partially equilibrated oscillators (e > e*) for which the equipartition theorem (40) holds out of 
equilibrium as in the model without disorder 32 ; on the other hand, there are non-equilibrated modes (e < e*) for 
which the equipartition theorem does not hold. This fact has been illustrated with the binary case of a system with 
two different values of e, i.e two different time-scales. 

In the next section we will define an effective temperature for each energy mode (each time-scale) in the system. 
Again we will distinguish between the partially equilibrated modes and the rest. 
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B. MODE CORRELATIONS AND RESPONSES 



We can define the correlation and the response functions for each energy mode. Let us define the mode-correlation 
C(e, t, s) as: 

C(e,t,s)=<±J26(e i -e)x i (s)x i (t)> , (41) 

i 

which is the correlation between the positions of the oscillators which have the same spring constant Ke. In this case, 
similar calculations as in 32 , show that the dynamical evolution of the mode-correlation function is given by: 

dC(e,t,s) ( dE t . A e 

We can compute the dynamic evolution for the mode-response after perturbing the Hamiltonian: 

H e = Ha — hM(e) , (43) 

where H is defined in (5) and the perturbation is coupled to the magnetization of that particular mode of energy e. 
We define the mode-magnetization as: 

M(e,t)=<±Y / 5(e i -e)x i {t)> ■ ( 44 ) 

i 

The mode-response function G(e, t, s) is defined as usual, 

G(e,t,s)=( S 4m t>s . (45) 



5h(a) 



Computing the dynamical equation for the mode-magnetization and deriving it respect to the external field h, we 
obtain the dynamical equation for the mode-response function: 

dG(e,t,s) (dE A . A e . f dE t . A g(e) r/ 

Note that the evolution of the mode-correlation and the mode-response are identical and differ only from a delta 
function term. Equations (42,46) can be solved exactly and allow us to make a deeper analysis of the non-equilibrium 
behavior of the system. The equilibrium behavior of these functions is analyzed in the Appendix B. 

In the general nonequilibrium regime the correlation function and the response function can be written (B2,B11): 

2 

C(e,t,s) = — w(e,s)exp[-(t- s)/T(e,t,s)} t > s, (47) 

G(e, t, s) = F(e ^ g(f) exp [-(t - s)/r(e, t, s)} t > s, (48) 
where F(e, s) is given in (B3) and we have defined a time-dependent relaxation time r(e, t, s) (B4): 

(t-s) ( 1 



r(e,t,s)='—^\ M ; | t>s. (49) 

Is 2Q 2 (x) a 



The relaxation time of a given oscillator is inversely proportional to the value of its spring constant e. This result 
indicates that each type of oscillator (each e) can be considered as a single variable (evolving according Langevin 
dynamics) with a time-dependent relaxation time r(e,t, s) given by the expression (49). Note that the distribution of 
disorder fully determines correlations and responses through the functions f{x) and Q2(x). 
The full correlation function is then given by, 



1 

C(t, s) =< — Xi{s)xi{t) >= J deC(e, t, s) 



(50) 
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(an analogous expression is valid for the response function) . The relaxation of the full correlation function depends 
strongly on the disorder distribution g{e). In equilibrium, the full correlation function can be written using (B7,B8,B9), 

^-'-n^--^))* ■ <5i) 

Defining the parameter r = ^ ^~^ K ^j this correlation (51) is given by the following Laplace transform, 

/•OO 

C oq (T)=/ h(e) exp (-er) de , (52) 
Jo 

where we have defined h(e) as: 

He) = ^ ■ (53) 

The equilibrium correlation function C eq (r) is well defined whenever the integral (52) converges, i.e. for g(e — > 0) ~ e a 
with a > 0. In this case, also the entropy constant in (7) is finite. When the integral (52) does not converge it means 
that the largest contribution to the full correlation function is determined by diffusive behavior of the e — > modes. 

A relevant feature of the DOM is that it can describe the heterogeneous behavior observed in glass formers (see 
the final discussion in Section VII). The decorrelation time in the DOM can be defined in terms of the integrated 
relaxation time 16 : 



J °°C^(T)dT tZ~ 2l /T%r eC 



Trclax - - ^ R - \j ^ R ^ CX P 



-2 



KA 2 e 1 
8 T 



(54) 



where we have used (20) (with e~ = ei) and the expression for the viscosity 7 (BIO). This expression shows that the 
relaxation time is thermally activated and that the main effect of the disorder on the relaxation time is through the 
activation barrier which is proportional to e. The other moments of the disorder distribution e~T, eZi only affect the 
prefactor of the relaxation time. 

In general, the distribution of disorder g(e) can be chosen in order to reproduce any given relaxation function by 
computing the inverse Laplace transform of the full correlation C cq (r). As an example we consider a few cases: 

• Binary distribution. Consider a binary system with one half of the oscillators having e = e\ and the other half 
having e = e 2 where e 2 > e\. In this case, the full correlation function in equilibrium (51) is simply the sum of 
two exponentials, 

C ° q (*- S )= 2^T 2 ( e2CXp ("^r (t_s) ) + £ i cx p(-^ a • (55) 

For long times, the full correlation (55) is dominated by the first exponential. For a general uniform discrete dis- 
tribution with n finite different values of e, the relaxation of the full correlation will be the sum of n exponentials 
of the type (B7). 

• Continuous distribution. Consider now the continuous distribution 17(e) = Ae a exp(— e) where A is the normal- 
ization factor. In this case, the full correlation function in equilibrium can be computed from (50,51) giving a 
power law decay: 

^'-'-^(h^j)" ■ < 56 » 

• Stretched exponential relaxation. Let us consider now the inverse problem of finding the disorder distribution 
for a correlation of the type: 

C cq (t- s) occxp(-AVt^) . (57) 
Taking the inverse Laplace transform, we get for the disorder distribution: 

1 -y 

g(e) oc e~2 e -47f7 . (58) 
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where A is a positive constant. Note that this distribution is non-normalizable since it decays too slowly in the 
e — > oo limit. This is not a problem since the large e limit only influences the short [t — s) behavior of the 
correlation function. Therefore by introducing a cutoff for g(e) in the large e limit we do not modify the stretching 
of the relaxation function. For the large-i generic stretching exponential behavior C oq (Y) <~ exp (— Ar b ^ we have, 

in the e — > limit, g(e) <~ exp(— B/e 1 ^) with B a positive constant. Taking b = 1/2 we recover the result (58) 
in the e — > limit. The exponential relaxation corresponds to the case 6 — > 1 and the g(e) has a gap. In general, 
the large t behavior of C cq (r) is determined by the e — > dependence of g(e). 

From this analysis we can see that the DOM is a very general model. The distribution of the disorder can be chosen 
in order to fit any experimental data just computing the inverse Laplace transform of the correlation. 



C. MODE EFFECTIVE TEMPERATURES 



Building up on the ideas by Cugliandolo, Kurchan and Peliti 9 we can define an effective temperature for each mode 
and analyze how the different degrees of freedom in the system relax. From (1) we can define an effective temperature 
for each mode: 



^-{<^)-W)\«* ,) + fm-s^) ■ (59 » 

The mode effective temperature T c g(e, s) only depends on the lowest or waiting time s but not on the largest time t. 
This is a typical feature of mean-field structural glass models 18 . 

As mentioned in Sec I, we can also define an effective temperature in the Fourier space (2) for each mode: 

T '«< eA »> = n |fetf • < 60 » 

where S(e : s) and x(e, Q, s) are denned as 16 : 

dtcoa(nt)C(e,t,s) , (61) 

dt\(e, t, s) exp (— iVtt) , (62) 

In equilibrium the power spectrum can be written as: 

/•oo 

5 oq (e,0)=/ drcos(OT)C eq (e,t-.s = r) , (63) 
where C cq (e, t, s) is given in (B7). The result is a Lorentzian function: 

S eq (e,0) = ^(e) ( 1— V J . ((,4) 
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The full power spectrum is expressed in terms of the full correlation function (50) and is given by the integral of the 
Lorentzian weighed by g(e): 

S^(Q) = ^ dr cos(ttT)C oq (r) = - T g(e) I i— - I de . (65) 

As for the full correlation function the disorder distribution g(e) can be chosen to fit the experimental power spectrum 
data by using the relation (65). In (62) we have used the susceptibility, which in terms of the response is given by, 

X (e,t,s)= f G(e,t,T)dr . (66) 

J s 
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Note that there is a factor 2 missing in (60) as compared to the standard definition of the FDT in frequency space 16 . 
This is due to the fact that we are assuming that the integral (62) extends from s (when the measure of the relaxation 
function starts) to oo rather than from — oo to oo as in the equilibrium case. The calculations are shown in Appendix 
D. With these definitions we obtain for the effective temperature in Fourier space (D6): 

T&(e,Sl,8)=T&(e,8) = -^w(e,8) . (67) 

Interestingly, the effective temperature in Fourier space does not depend on the frequency and only depends on the 
waiting time s and the mode e considered. We can observe that the result in Fourier space (67) is equal to the first 
term of (59) in the time domain so the subleading corrections introduced by the second term in (59) are absent in 
Tf s (e,s). _ 

Let us note that it is also possible to define a "full" effective temperature T e ff(s) (either in the time or frequency 
domains) from the full correlation function (50) and the analogous expression for the response function. From the 
analysis carried out in Section IV A we have that e* — > asymptotically. The full effective temperature is then 
determined by the contribution to the full correlation and full response from the modes in region I, e > e*. In Fig. 5 
we can see the time-evolution (in Monte Carlo sweeps) of the normalized full effective temperature in time domain 
for the system considered in Fig. 3, with 100 types of oscillators with spring constants uniformly distributed in the 
range e = 0.01 — 1. As time goes on, more and more oscillators partially equilibrate and the full effective temperature 
approaches asymptotically to 0(s) (27). 



0.93 




FIG. 5. Time evolution (measured in Monte Carlo steps) of the full effective temperature (normalized by its asymptotic 
value (27)) for the same system considered in Fig. 3. 



In the glass literature it is sometimes useful to construct the so-called fluctuation-dissipation (FD) plots 40 . Generally, 
a FD plot is a plot of the integrated-response versus the correlation function. In the DOM, the full correlation (50) 
at equal times C(s,s) depends on s. Therefore, for a meaningful FD plot we should consider the normalized full 
correlation C(t, s) = C(t, s)/C(s, s) and the full integrated-response given by: 

/t />oo 
dr \ deG(e,t,r) . (68) 
Jo 

The FD plot consists on representing \ as a function of C for a fixed waiting time s and varying t. At equilibrium 
this plot gives a straight line with slope — 1/T where T is the temperature of the bath. In the DOM, for asymptotic 
waiting times (s — > oo) we find, for a fixed s, 

X(5) = ^(l-e) • (69) 
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Therefore, for a given asymptotic waiting time s the FD plot gives a straight line with negative slope equal to the 
inverse of the asymptotic effective temperature (27). In Fig. 6 we represent the resulting asymptotic FD plot which 
has also been obtained before in the OM model 36 as a particular case. Straight FD plots are characteristic of the 
one-step behavior of structural glasses. In Fig. 6 we can also see the absence of an equilibrium relaxation (usually 
called /3-relaxation) for small responses. This unrealistic feature is due to the mean-field nature of the model. 



S > oo 







c 



FIG. 6. Asymptotic FD plot: the full integrated response function normalized by the asymptotic effective temperature (27) 
as a function of the normalized full correlation (69) for s fixed. 



Next we consider in detail the effective temperature for the different types of modes. As we will see all modes in 
region I give the same mode effective temperature value which is equal to the value of the full effective temperature. 



1. Region I of partially equilibrated modes: e > e* 



For the partially equilibrated modes, i.e e > e*, we just keep the first term in (59) (the derivative of the w(e, s) goes 
to zero faster than w(e,s) itself). Therefore, 



T cff (e, s) ~ ^y^(e, s) with ew{e, s) = ^Q 2 (s) , (70) 



and the effective temperature reads: 



T cff (e, S )~ |q 2 («) = 6(s) , (71) 



which coincides with the expression derived within the adiabatic approximation (27) (up to subleading logarithmic 
corrections of order l/log(t)). 

For the effective temperature in Fourier space we get the same result but without logarithmic corrections: 

Tf s (e,s) = -^ )W (e,s) = Q(s) . (72) 

All modes which are partially equilibrated share the same effective temperature that we define as 0(s). As the time s 
goes on, more and more oscillators achieve this value Q(s) and, in the asymptotic limit, all oscillators will thcrmalizc 
to this effective temperature. 
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FIG. 7. (Color online) The dependence of the effective temperature (in time domain) as a function of the value of e for the 
same system as in Fig. 3 with 100 types of oscillators with e- values equally spaced and uniformly distributed between 0.01 and 
1.00. Data is shown for different waiting times at zero temperature. In the inset we show the scaling of the effective temperature 
with the threshold energy e*. 

In Fig. 7 we can see the dependence of the effective temperature (59) on the energy mode at a given time. We plot 
the effective temperature normalized by 0(s) for different values of s at zero temperature. As we increase the waiting 
time s those oscillators with higher value of e (shorter relaxation time according to (B8)) tend to achieve an effective 
temperature T c g(e,s) equal to 0(s). Whereas those oscillators with lower values of e have an effective temperature 
greater than 0(s). From the figure we note that the ratio T e g(e, s)/0(s) for e > e* is smaller than 1 due to the 
logarithmic correction coming from the second term in the r.h.s of (59) whereas such corrections are not present in 
T^ s (e, s). Therefore, as the waiting time increases, lower energy modes progressively reach partial equilibrium. In the 
inset of Fig. 7 we plot the normalized effective temperature as a function of e/e* showing good collapse. The shape 
of Fig. 7 is expected to be quite general. The same qualitative picture was found for a Gaussian distribution for the 
disorder or an exponential one (data not shown). 

Modes around the threshold, e ~ e*, are not yet partially equilibrated. By plotting dw ^ s ^ as a function of e we can 
identify the width Ae of the region separating the partially equilibrated modes from the rest. This is shown in Fig. 8 
for the same system considered in Fig. 7. There, we have multiplied dW Q € s ' s ^ by a factor log S ( s ) which is the inverse of 

the leading order of the decay of dW g t s ' s ^ computed in the adiabatic approximation. In the inset we can clearly see 
that this quantity again scales as a function of . From such scaling we note that the width Ae around the threshold 
scales as e* and therefore, as e* decays to 0, the region around the threshold becomes progressively thinner relative 
to the support of g(e). This allows us to define the threshold boundary as the asymptotically thin sector of e- values 
that separates region I from region II. 
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2. Region II of diffusive modes: e < e* 
Using (30), the expression for the effective temperature (59) can also be written as: 

T off ( e , s )=26( S )-^M . (73 ) 
.9( e ) 

As a consequence those degrees of freedom with larger time-scales are hotter than partially equilibrated ones. This 
result can be checked in Fig. 7 where the oscillators in the region e < e* have an effective temperature higher that 

For finite times, in the limit e — > 0, the second term in the rhs of (73) vanishes 

T cS {e -> 0) = 20(a) , (74) 

and the effective temperature is just twice the asymptotic effective temperature (the effective temperature for the 
partially equilibrated modes). Note that this value corresponds to the effective temperature of a diffusive process , 
whose value is twice the temperature as given by the Einstein relation D = 2T where D is the diffusion constant 
associated to the diffusive variable x, (Ax 2 (t)) ~ Dt. This result can be intuitively understood. Modes with e small 
do not contribute much to the total energy (5), so their motion is completely uncorrelated from the rest of modes and 
become diffusive. 

In order to illustrate this result we show in Fig. 9 the effective temperature for a binary system. In this case, 20% 
of the oscillators have e = 0.1 while 80% have e = 1. We can see how the oscillators with larger e are partially 
equilibrated to the effective temperature Q(s) (71) while the oscillators with smaller e have an effective temperature 
approximately equal to 2Q (s) as given by (73). 
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FIG. 9. We plot the effective temperature for a binary system: 20% of oscillators have e = 0.1 (triangles) while the rest have 
e = 1 (circles). We also plot the time-evolution, measured in Monte Carlo steps, of O(s) and 20(s) (71) (full lines) for the two 
modes. 



V. FLUCTUATION RELATIONS AND THE EFFECTIVE TEMPERATURE 



The effective temperatures for the partially equilibrated modes can be obtained by using a recently proposed 
fluctuation relation . In it has been suggested that effective temperatures can be derived by assuming microscopic 
reversibility between surfaces of constant energy. This relation has been also proposed in the framework of mean-field 
spin glass models 44 and applied to generalized oscillator models 36 . The effective temperature O(s) can be computed 
from the ratio between the probabilities of having a change in energy equal to SE and (-8E). From (11), 



PjSE) 
P(-SE) 



= exp SE 



2Q 2 (s) 



= exp 



SE 



(75) 



where we have used (11). This relation between forward and reverse transitions only holds at the level of energies rather 
than configurations. It should not be confused with the detailed balance or reversibility property (4) of microscopic 
dynamics which obviously holds at the level of configurations in and out of equilibrium (4). Equation (75) defines 
a nonequilibrium temperature that is equal to the effective temperature obtained from the fluctuation-dissipation 
relation for the partially equilibrated modes (67,71). 

We can also extend this argument and compute the probability P s {Q e ) of having an instantaneous change in the 
energy Q t of a given mode e at a given time s, 



Q, 



(76) 



The probability P s (SQ e ) is computed in (A9). It is given by: 

1 



I'^Q.) = - 77 , 

,JA-kK A 2 e 2 w(e, s) 

and the fluctuation relation for each mode is given by: 



exp 



P s (SQ t ) 
Ps(SQ e ) 



exp 



9(e) " 
2ew(e, s) 



(We 



4KA 2 ew(e, s) 



= exp 



(77) 



(78) 
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Only for partially equilibrated modes, we can use (71) and get, 



Ps(SQ e ) 
Ps(SQ e ) 



exp 



SQ e ' 
G(s)_ 



(79) 



This relation emphasizes the fact that only for the partially equilibrated modes a unique effective temperature can 
be properly defined. Otherwise, for non-equilibrated modes 2£ ™^'^ ^ 0(s) 

Recent work 36 ' has shown that this effective temperature can also be obtained from microcanonical relations 
relating observable changes. In the present case, we consider the mode- magnetization (44). The probability of having 
a change 5M t in the mode- magnetization is computed in the Appendix C (C6). The rate of having a positive and a 
negative change gives: 



W h (SM e ) W°(SM e ) / e \ W°(SM e ) ( SM e h\ 

J ■ exp 5M e h—— = — ^ exp — . (80) 



W h {-5M e ) W°{-5M e ) ^\ f - 2Q 2 (s)J W°{-5M e ) ^\0(s) 

At difference with (78) now the effective temperature O(s) describes mode- magnetization fluctuations for all values 
of e. This is not unexpected due to the linear dependence of the observable M e on the configurational variables (44) 
and the Gaussian character of the distribution (CI) which determine that (C6) depends on Q 2 (i.e. on 0(s)) for all 
values of e. Although not proven we suspect that this general dependence might not hold in general, e.g. for generic 
non harmonic potentials 36 . 



VI. SUMMARY OF RESULTS 



In this paper we have introduced an exactly solvable disordered model for glassy dynamics, the disordered oscillator 
model (DOM). The model is the disordered version of the oscillator model (OM) previously introduced by one of us 32 
as a simple model with entropic barriers for aging dynamics. The DOM consists of a collection of non-interacting 
oscillators or modes with different spring constants. The model has trivial statics but is kinetically constrained and has 
cooperative dynamics displaying glassy behavior at low temperatures. In the DOM the relaxation time associated to 
each mode is inversely proportional to the effective spring constant. This allows us to study the relaxational properties 
(both in equilibrium and in the aging regime) associated to each mode and the value of the effective temperature (1) 
for the different degrees of freedom. 

Dynamics in the DOM can be solved writing down a hierarchy of coupled first-order differential equations (19). 
These set of equations can be closed by introducing a generating function and using the method of characteristics or 
the adiabatic approximation at low temperatures. The latter approach is much simpler and allows us to define an 
effective temperature O(s) (27) in the aging regime (s being the waiting time). The relaxation of one-time quantities 
appears similar to the OM, just replacing the spring constant in the OM, K, by the average spring constant Ke in 
the DOM. 

To investigate how different energy modes relax to equilibrium we introduce the mode-energy density (28), the 
mode- magnetization (44) and the mode-correlations and mode-responses (41,45). From the dynamical equations 
equations (30) we can prove the existence of time-dependent threshold located at e*(s) separating two regions (I and 
II) characterized by the sign of the time derivative of the mode-energy density (35,36,37). Only the modes in region 
I satisfy the equipartition relation at the effective temperature O(s) (40). The threshold value e*(s) decays with the 
age of the system with more modes entering region I as the system approaches equilibrium. In Fig. 10 we show an 
schematic picture of the different mode regions in the model. 

Furthermore, for each mode we can define an equilibrium relaxation time (B8) in terms of a temperature dependent 
viscosity that diverges when T — > (BIO). Full correlation functions diverge if too many diffusive modes are present, 
i.e. if g(e) vanishes too slowly in the limit e — > 0. For any equilibrium correlation function C oq (r) there exists a 
distribution g(e) that in the limit e — > reproduces the long-time decay of C cq (r). 

Mode-dependent effective temperatures can be also defined (59) from mode-correlations and mode-responses. They 
agree with the previous scenario. In region I modes are mutually equilibrated at the temperature 0(s) whereas 
in region II modes are out-of-equilibrium. The latter show diffusive behavior and are characterized by an effective 
temperature larger than 0(s) and equal to 20(s) for e — > 0. This last result is characteristic of purely diffusive behavior 
(Xoo = 1/2, see 16 ' 18 ' 41 ). The effective temperature can also be computed in frequency space (2) corroborating the 
results obtained in the time-domain. It is noteworthy to mention that no frequency dependence is found (67). This 
result is in contradiction to what is observed in experiments and might be consequence of the simplicity of the model 
(see the discussion in the next section). 
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Finally, we have corroborated the scenario of mode-dependent effective temperatures using the fluctuation relations 
recently proposed in the framework of aging systems 36,42,43 . Interestingly, for both partially equilibrated modes 
(region I) and diffusive modes (region II) the fluctuation relation for the mode-magnetization gives the unique value 
O(s). It is unclear to us whether this result is of general validity or restricted to the DOM. It would be interesting to 
put under scrutiny other types of models. 



REGION II Ae REGION I 




NON— EQUILIBRATED PARTIALLY 

EQUILIBRATED 



FIG. 10. (Color online) Schematic picture showing the two regions of modes. Region I (e > e*) consists of partially 
equilibrated modes at the temperature O(s). In region II (e < e*) there are the non-equilibrated or diffusive modes. In between 
the two regions there is a boundary region of width Ae around the threshold e* where modes start to equilibrate. As time goes 
on, the threshold value e* decreases with the width Ae being proportional to e*. For sake of clarity the size of the regions in 
the picture is not in scale (region II and Ae should be tiny as compared to region I). 



VII. CONCLUDING REMARKS 



In which ways can the DOM improve our current understanding of the glassy state? A lacking piece in the glass 
transition puzzle is to understand whether it is possible to identify effective parameters (such as local temperatures or 
stresses) that reflect the nonequilibrium relaxation toward equilibrium of macroscopic or bulk properties (e.g. enthalpy, 
specific heat, volume, pressure, viscosity, thermal conductivity,...). However, if dynamics and transport phenomena 
are strongly heterogeneous then it is hard to understand how a description of the macroscopic time evolution of bulk 
quantities can be achieved from a few number of parameters as proposed by several phenomonological approaches. 

It would be very interesting to get experimental evidence that a unique nonequilibrium or effective temperature 
can be identified in glass formers in the aging state. Such temperature would constitute a candidate parameter to 
describe many physical nonequilibrium properties in the aging state. This point of view has been advocated a few 
years ago by Nieuwenhuizen 45 . The existence of a unique effective temperature appears to be an inherent property of 
mean-field systems. The DOM is yet another example of a mean-field system with a well defined nonequilibrium or 
effective temperature (that we have denoted by 0(s)). Nevertheless, the current experimental evidence points toward 
the fact that dynamics is strongly heterogeneous and a unique nonequilibrium temperature might not exist. Rather, 
the effective temperature appears more as a local quantity that changes from region to region inside the glass and 
which takes a value that depends directly on the degree of mobility of the region. 

The DOM is a mean-field model with cooperative glassy dynamics. Alone it cannot explain the heterogenous 
behavior observed in glasses because it does not incorporate the fact that different modes intermittently switch for 
mobile to frozen behavior. The DOM describes the dynamics of an isolated cooperatively rearranging region (CRR) 
maybe of a few nanometers of spatial extension that includes a few hundreds of atoms or molecules. The different 
energies e of the vibrational modes in each CRR are determined by its volume (the number of atoms it contains) and 
its surface through possible interactions with nearby CRRs. The glass could then be viewed as a mosaic made out 
of many CRRs each one described by a DOM with a given energy levels distribution g(e). Such distribution strongly 
fluctuates among CRRs. In the scenario where CRRs intermittently appear and dissapear at a given rate, each time a 
new CRR is produced it carries a new energy level distribution g(e). Consequently, at a given time each region shows 
a different e. According to (54), fluctuations in ~e will lead to a great disparity of relaxation times among CRRs and to 
the presence of slow and fast regions. Also, the effective temperature will be inversely proportional to e~ (71) leading 
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to an effective-temperature age-dependent field 0(e,s) as well as a ^-dependent effective temperature, T c ^(f2, e, s). 
The latter dependence of the effective temperature on the probing frequency would agree with the experimental 
results 11 ' 31 . A schematic picture of this scenario is shown in Fig. 11. This microscopic scenario is similar to that 
proposed by Xia and Wolyness 22 from the random first order transition theory for mean-field spin glasses 15 ' 46 applied 
to the glass transition. A related scenario in terms of aggregating activated regions has been also described by Crisanti 
and Ritort 47 . A more elaborated and compelling description of the mosaic picture has been recently proposed by 
Bouchaud and Biroli 23 . 

Merging of random energy levels and cooperative dynamics in spatially localized regions (CRRs) appears a promising 
route to provide a phcnomcnological description of the glassy state. It remains to be seen whether such view is capable 
to provide new predictions that can be tested in the laboratory. This maybe difficult until we get a more clear picture 
and better understanding about the nature of the heterogeneities uncovered by the experimentalists. In turn, this 
might provide a definitive experimental evidence in favour of the intuitive notion of CRRs put forward by Adam- 
Gibbs-Di Marzio many years ago. 




FIG. 11. Schematic picture of a piece of glass made out of different CRRs, each one characterized by a value of e (we show 
9 labelled regions in the picture). In such a mesoscopic description regions are maybe of a few nanometers of spatial extension 
and typical activation energies KA 2 ~e are of the order of a few ksT. lvalues fluctuate from region to region. According to (54) 
dynamics is highly heterogenous with relaxation times in local regions differing by many orders of magnitude. The effective 
temperature (71) is also inversely proportional to ? leading to an effective-temperature field. 



APPENDIX A: DYNAMICAL EQUATIONS FOR Q K 

To solve the hierarchy of equations for Qk we start by computing the evolution for the quantity Q2W (12) which 
has appeared in the dynamic evolution for the energy (14). In order to do that we compute the joint probability 
distribution of having a change of SE in the energy (5) and 8Q2 in Q 2 (12) for an elementary move (8), 



P(5E,SQ 2 )= [ 6[5E-KJ2< 



TiXi r 

+ 2iV, / 



This probability factorizes and can be written as follows: 



P { *,y) = PW^^(dM^ , (A2) 
Ji = K^fl + {x _ ac) Qs > j 2=K A 2 (^+Q,), (A3) 



2 v 'Q2 ' \Q2 
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where x = SE, y = SQ 2 and P(x) is given in (11). The other quantities can be defined through the general expressions: 

i i 

From the joint probability we can derive the equation for Q2(t) using (4): 



dQ 2 
dt 



dx dyyP(x,y) + dxe~ &x / dyyP(x,y) . (A5) 

-oo J —oo JO J —oo 



The integral over y is just the quantity J\ and does not depend on J 2 . Expressed in terms of the dynamical evolution 

dE 
at 



of the energy W- and the acceptance A we find: 



dQ 2 (KA*e 2 a cQAA{t)+ (Qsy_E (Ag) 



dt V 2 Qi ) \Qi) dt 

As we expect, the evolution of Q 2 depends also on the new quantity Q3. Therefore we have to generalize these 
computations in order to obtain the whole set of coupled equations. The result is: 

These equations are, for a given probability distribution of the disorder, the complete solution for the dynamical 
one-time quantities. 

At zero temperature only negative changes in energy are accepted. The equations turn out to be a bit simpler: 



d Qk _ K6?t k _ Q^lkVQ* ^ _ 



dt 2 v ' Q 2 V 7T 4fcA 2 Q 2 

We now outline the computations for the mode-energy density. As for the quantities Qk we must first compute the 
joint probability of having a change of SE in the energy (5) and Sw in w(e) (28), the result is: 

^--^t^M^^F) • (A9) 

where now x = SE, y — Sw and P(x) is given in (11). The quantities Si(e) and 62(e) are given by, 

S 2 (e) = KA 2 w(e) - i^A 2 ^ , (A10) 

Q2 

ew(e) . . KA 2 g(e) . , , 

Si(e) = -A±( x -a c ) + f-^ . (All) 



From the equation, 



dw(e, t) 



/0 poo roc />oo 

dx dyyP(x,y)+ dxe~ fix dyyP{x,y) , (A12) 
-00 J — 00 Jo J— 00 



dt 
we get: 

dw(e,t) few(e,t)\ dE (K A 2 g(e) a c ew(e,t) 



dt \ Q 2 (t) J dt V 2 Qi{t) 



A(t) . (A13) 
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APPENDIX B: EQUILIBRIUM BEHAVIOR OF THE CORRELATIONS AND RESPONSES 



In this appendix we will solve the dynamical equations for the mode-correlation and the mode-response functions 
and analyze their equilibrium behavior. First of all, let us define: 

f(t) = -(j^-a c A(tfj . (Bl) 
In terms of this function the solution for the mode-correlation (42) is given by: 

C{e,t,s) = C{e,s,s)e^[-e J^^-^dx^j . (B2) 
Where the initial condition for the mode-correlation is: C(e, s, s) — ^w(e, s). In a more compact way we can write: 



C(e,t, s) ~ -^w(e, s) exp ^— J F(e,x)dx^j with F( 



<- c)=*4tL ■ (B3) 

2Q 2 (x) 



In general we can define a time-dependent relaxation time, 



r(e,t,s)= { -L-A\ - t - - 1 (B4, 



Js 2Q 2 (x) aX . 



and the correlation function can be written as: 



C(e,t,s) = ^w(e,s)cxp[-(t-s)/T(e,t,s)} . (B5) 

In equilibrium is given by (21) and / eq (Bl) is given by: 

/ oq = a c crfc(a oq ) , (B6) 
where a cq is given in (15). Therefore, the mode-correlation in equilibrium has an exponential behavior: 

2?/) cq 

C oq (e,t, S ) = -^ r exp[-(t- S )/r oq ( e )] , (B7) 
where w eq is given in (31) and the characteristic equilibration time of the energy mode e is: 

- h (a%§^)) ■ < B8 > 

As we have mentioned before the characteristic relaxation time of all the oscillators belonging with energy mode e 
is inversely proportional to the value of e. Therefore, the higher is the value of e the shorter is the characteristic 

time-scale of the oscillator. Therefore we can define the friction coefficient 7 using the relation r cq (e) = for an 
overdamped harmonic oscillator in thermal equilibrium: 

7 = (A 2 erfc(a cq )) ' (B9) 

which, for low enough temperatures behaves, 

7=\/^r ex P -5- ^ • ( B1 °) 
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Therefore, the divergence of the relaxation time is of an activated type just as in the model without disorder. 
We can also derive the equilibrium relaxation time from the response function: 

G(e, t, s) = giflfMf) exp (- J* F(e, x)dx^j 6(t - s) , (Bll) 

with F(e,x) given in (B3). In equilibrium the mode-response becomes: 

GCq(e ' *' s) = 7$$jK exp [_( * ~ s)/r ° q(e)] 9{t s) ■ (B12) 

Now we can check that in equilibrium FDT is verified : 

/ dC(e,t,s) \ „ / n 

I \ 2«,(e) T (B13) 



\G(e,t,s)J g(e) 

APPENDIX C: MODE-MAGNETIZATION RATES 

In this section we derive the main results for the mode-magnetization rates. The joint probability of having a 
change in the mode magnetization (44) 6M € and in the energy 5E can be computed as in the original model 34 ' 43 . The 
result is: 



47r^5 2 (e)5 4 (e) V 4S 2 (e) / " ^ 4S 4 (e) 
where we have defined the quantities: 



5 2 (e) = KA 2 Q 2 ~heKA 2 M + <D{h 2 ) , (C2) 
5 3 (e)=£MKA 2 -/iA 2 9 (e) , (C3) 

S 4 (e) = ^|^+0(/> 2 ) ■ (C4) 
At T = the probability that an attempted change 5M e is accepted is given by: 

,o 

W h {6M e )= P h {5E,5M e )d(5E) . (C5) 

»/ — oo 

Assuming that the value of the magnetization is zero (condition of neutrality ' ) we find: 



W h (5M e ) = W (5M e )exp \bM f h-{— 



(C6) 



where 



W°(SM.) = exp U™£) f exp ( ^f-f ) ^E) , (C7) 

4ny/S2(e)SA(e) V 4S 4 (e) J J _ x y \ 45 2 (e) 1 



and therefore, 

W h (tfM £ ) 



W h (-<JM £ ) ^ V 2Q 



expj^/^-L-) (OS) 
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APPENDIX D: THE EFFECTIVE TEMPERATURE IN FOURIER SPACE 



The effective temperature in Fourier space is denned as (60): 

^•"•*> = n TOi ■ < m » 

fn our case, the real part of the power spectrum is given by, 

2 f°° 

S(e,n,s) = -w(e,s) J dt cos(Ot) exp[-(t - s)/r(e, t, s)] , (D2) 
where we have used (47). The susceptibility is, 

X (e, t,s)=J G(e, t, r)dr = ^ (1 - exp[-(i - s)/r(e, t, a)]) , (D3) 
where we have used (48). The complex part of the response is defined as (62): 

dt X (e,t,s)e- mt =x'(e,n,s)-ix"(e,n,s) . (D4) 

Taking the imaginary part we get, 

X"(e,^,s) = ^^ dt cos fit exp[-(t - s)/r(e,t,s)] . (D5) 

The final result is: 

Tf s (e,ft,s) = ^w(e,s) . (D6) 

Note that the effective temperature in Fourier space does not depend on the frequency Q. This result is consequence 
of the identical functional decay of the correlation and response functions (B3,B11) which leads to nearly identical 
expressions for the power spectrum (D2) and the frequency-dependent susceptibility (62). In this respect this result 
is not generic of mean-field models which instead do show a dependence of the effective temperature T g(e, fi, s) on 
the probing frequency fl. 
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